%The scripts below are written in August 2008 to understand 
%1) The relation between CEJ as identified by negative EEJ values from
%CHAMP and the IEE data - done
%2) To verify the CEJ idetification with HUA magnetic data - done 
%3) Finally to produce a average CEJ signature from CHAMp and IEF data -
%ran into tourble. 19-aug-2008 - cmaaannoojjj

%%
data=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_scalar.txt');
%data=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_vector.txt');
%data=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\SAC-C_scalar.txt');
%data=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\SACOERSTEDCHAMP_Scalar.txt');
load c:\manoj\projects\plasma\HUA_PIU1.mat HP HUAPIU HUA_PIU_fday;
load c:\manoj\projects\ace\OMNI_ELEC ace_all ace_gsm_bz; % JUST ACE IEF
load c:\manoj\projects\plasma\HUA_PIU1.mat magnetic_frac;
% load('C:\Manoj\projects\eej\YAPMNBIK.mat');
% HP(1,:) = yapmnbik;
% HP(2,:) = yapdate;
% HUAPIU=YAPMNBIK;
% HUA_PIU_fday = unique(floor(yapdate));
ltt = champ_lt_new(data(:,1),data(:,2));


%%
L = data(:,2) > -90 & data(:,2) < -60 & lt > 9 & lt < 18; %HUA
%L = data(:,2) > 125 & data(:,2) < 155 & lt > 9 & lt < 18; %YAP

mydata = data(L,:);

% kp=1;
% for i = 1:length(mydata),
% [trash ind] = min(abs(HP(2,:) - mydata(i,1)));
% dummy1(kp) = ind;
% tdifer(kp) = trash;
% kp=kp+1;
% end;
% L = tdifer < 0.0014;
% huadata = HP(1,dummy1(L));


win_dow= 10/(60*24);%10 minutes before and after

kp=1;
for i = 1:length(mydata),
a=find(HP(2,:)<(mydata(i,1)+win_dow) & HP(2,:)>(mydata(i,1)-win_dow));
if any(a)
dummy(kp) = a(1);
huadata(kp) = nanmean(HP(1,a));
index_h(kp) = a(1);
index_i(kp) = i;
kp=kp+1;
end;
end;



seldata = mydata(index_i,5);
kpdata = mydata(index_i,10);
champtime = mydata(index_i,1);
huatime = HP(2,index_h);

LL = seldata < 0;

seldata = seldata(LL);
champtime = champtime(LL);
huatime = huatime(LL);
kpdata = kpdata(LL);



for i = 1:length(huatime),
    ia(i)= find(HUA_PIU_fday==floor(huatime(i)));
end;

%[c,ia,ib] = intersect(HUA_PIU_fday,floor(huatime));

meanhuadata = nanmean(HUAPIU);


for i = 1:length(huatime),
        L = ace_all(:,1) >= HUA_PIU_fday(ia(i))& ace_all(:,1) <= HUA_PIU_fday(ia(i))+1;
        plot(floor(huatime(i))+magnetic_frac, meanhuadata,'r');
        hold on;
        plot(ace_all(L,1),ace_all(L,2)*10,'k');
        plot(floor(huatime(i))+magnetic_frac,HUAPIU(ia(i),:));
        plot(champtime(i),0,'r*');
        title([datestr(champtime(i)+datenum(2000,1,1),21) sprintf(' kp=%5.2f champ cd = %5.3f ',kpdata(i), seldata(i)) ]);
        locp = find(data(:,1)==champtime(i));%location of champ pass in data
        plot(data(locp-5:locp+5,1),data(locp-5:locp+5,5)*1000,'k-*','LineWidth',2);
        hold off;
        pause;
end;









%%



champ=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_scalar.txt');
oersted = load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\Oersted_scalar.txt');
sacc = load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\SAC-C_scalar.txt');

load c:\manoj\projects\plasma\HUA_PIU1.mat HP HUAPIU HUA_PIU_fday;
load c:\manoj\projects\ace\OMNI_ELEC ace_all ace_gsm_bz; % JUST ACE IEF
load c:\manoj\projects\plasma\HUA_PIU1.mat magnetic_frac;
meanhuadata = nanmean(HUAPIU);


champlt = champ_lt_new(champ(:,1),champ(:,2));
oerstedlt =champ_lt_new(oersted(:,1),oersted(:,2));
sacclt = champ_lt_new(sacc(:,1),sacc(:,2));



L = champ(:,2) > -90 & champ(:,2) < -60 & champlt > 9 & champlt < 18; %HUA
champdata = champ(L,:);
L = oersted(:,2) > -90 & oersted(:,2) < -60 & oerstedlt > 9 & oerstedlt < 18; %HUA
oersteddata = oersted(L,:);
L = sacc(:,2) > -90 & sacc(:,2) < -60 & sacclt > 9 & sacclt < 18; %HUA
saccdata = sacc(L,:);

%%
number_of_days = 1;

for i = 1:length(HUA_PIU_fday),
    
    
      
    
    champ_i = find(HUA_PIU_fday(i)==floor(champdata(:,1)));
    oersted_i = find(HUA_PIU_fday(i)==floor(oersteddata(:,1)));
    sacc_i = find(HUA_PIU_fday(i)==floor(saccdata(:,1)));
    
    cej_check = any(champdata(champ_i,5)<0) | any(oersteddata(oersted_i,5)<0)| any(saccdata(sacc_i,5)<0);
    kpvalue = nanmean([mean(champdata(champ_i,10)) mean(oersteddata(oersted_i,10)) mean(saccdata(sacc_i,10))]);
%    if any([champ_i;oersted_i;sacc_i]),
    if cej_check && kpvalue <=2.1,
    
     win_dow =  rem(number_of_days,8);
    if win_dow==0,
        win_dow=8;
    end;      
    subplot(4,2,win_dow);      
        L = ace_all(:,1) >= HUA_PIU_fday(i)& ace_all(:,1) <= HUA_PIU_fday(i)+1;
         plot(HUA_PIU_fday(i)+magnetic_frac, meanhuadata,'r'); hold on;
         plot(ace_all(L,1),ace_all(L,2)*10,'k'); 
         plot(HUA_PIU_fday(i)+magnetic_frac,HUAPIU(i,:));
         
         if any(champ_i),
             locp = floor(mean(champ_i));
             locp = find(champ(:,1)==champdata(locp,1));
             
         plot(champ(locp-5:locp+5,1),champ(locp-5:locp+5,5)*1000,'k.-','LineWidth',1);
         end;
             if any(oersted_i),
             locp = floor(mean(oersted_i));
             locp = find(oersted(:,1)==oersteddata(locp,1));
               plot(oersted(locp-5:locp+5,1),oersted(locp-5:locp+5,5)*1000,'g.-','LineWidth',1);
             end;
             if any(sacc_i),
                 plot(saccdata(sacc_i,1),saccdata(sacc_i,5)*1000,'r.-','LineWidth',1);
             end;
             date_string = datestr(HUA_PIU_fday(i)+datenum(2000,1,1),29);
            title([date_string sprintf(' kp= %3.1f ',kpvalue)]);
            axis([HUA_PIU_fday(i) HUA_PIU_fday(i)+1 -inf inf]);
            %datetick('x','HH');
            hold off;
          
          number_of_days = number_of_days+1;
          
          if win_dow == 8 | i== length(HUA_PIU_fday),
            % h=legend('Observed JULIA EF (mV/m)','Predicted JULIA EF (mV/m)','ACE Ey (mV/m)/15','Location','Best');
            % set(h,'FontSize',8);
            set(gcf,'PaperPosition',[0.1,0.1,9,11.3]);
             print(gcf,'-append','-dpsc','-r600','-painters',['C:\Manoj\projects\cej\plots\HUA_CEJ']);
           % saveas(gcf,['C:\Manoj\projects\cej\plots\' date_string],'pdf');
            % display(['c:\manoj\projects\ace\plots\' date_string]);
            clf;
            
          end;
         
end;    
end;

  set(gcf,'PaperPosition',[0.1,0.1,9,11.3]);
  print(gcf,'-append','-dpsc','-r600','-painters',['C:\Manoj\projects\cej\plots\HUA_CEJ']);
        
%%
        lp = 1;
    
%     
%     for i = 1:length(HUA_PIU_fday),
%     
%     champ_i = find(HUA_PIU_fday(i)==floor(champdata(:,1)));
%     oersted_i = find(HUA_PIU_fday(i)==floor(oersteddata(:,1)));
%     sacc_i = find(HUA_PIU_fday(i)==floor(saccdata(:,1)));
%     
%     cej_check = any(champdata(champ_i,5)<0) | any(oersteddata(oersted_i,5)<0)| any(saccdata(sacc_i,5)<0);
%     
% %    if any([champ_i;oersted_i;sacc_i]),
%     if cej_check,
%         my_index(lp) = i;
%         lp=lp+1;
%     end;
%     end;
%     
%     
clear;
champ=load('C:\Manoj\projects\cej\stefan\EEJ_parameters_V3.0\EEJ_parameters_V3.0\CHAMP_vector.txt');
load c:\manoj\projects\ace\OMNI_ELEC ace_all ace_gsm_bz; % JUST ACE IEF
champlt = champ_lt_new(champ(:,1),champ(:,2));


L = champlt >9 & champlt < 18 & champ(:,10) <=2.1 & champ(:,5)> 0;
% L = champlt >9 & champlt < 18 & champ(:,10) <=2.1;% & champ(:,5) < 0;
%L = champ(:,10) < 4.1;
%L = champlt >9 & champlt < 14 & champ(:,5) < 0;

%L = champ(:,5) < 0;% & champ(:,10) <=2.1;

clt = champlt(L);
cut = champ(L,1);
cdt = champ(L,5);

clear cdata ctime ace_data_m;
numd = 1;
cpr=0;

for ii = 1:length(cut),
lodp = find(champ(:,1)==cut(ii));
if (cut(ii)>cpr&lodp>5),
     if mean(diff(champ(lodp-5:lodp+5,1))) > 0.06,
cdata(numd,:) = champ(lodp-5:lodp+5,5);
ctime(numd,:) = champ(lodp-5:lodp+5,1);
numd=numd+1;
%   else,
% cdata(numd,:) = champ(lodp-10:2:lodp+10,5);
% ctime(numd,:) = champ(lodp-10:2:lodp+10,1);
% numd=numd+1;
      end;   
cpr = champ(lodp+5,1);
end;
end;

 nn=1;  
for i = 1:length(ctime),
[trash ind] = min(abs(ace_all(:,1) - ctime(i,6)));
if abs(trash) < 0.0042, %(0.0042 ~= 6/(24*60)
ace_data_m(nn,:) = ace_all(ind-288:ind+288,2);
nn=nn+1;
end;
end;



 nn=1;  
for i = 1:length(cut),
[trash ind] = min(abs(ace_all(:,1) - cut(i)));
if abs(trash) < 0.0042, %(0.0042 ~= 6/(24*60)
ace_data_m(nn,:) = ace_all(ind-288:ind+288,2);
nn=nn+1;
end;
end;


% plot(mean(cdata),'k*')
% 
% figure(2);
% 
% for i = 1:577,
% plot(ctime_n(i,:)-ctime_n(i,6)+(cdata_n(i,6)/slope(i)),cdata_n(i,:),'k.');
% hold on;
% end;
% 

%%
clear cej_data
clear cej_array
clear cdata ctime
incrlat  = 30;
st=-180;
en=150;
nbin=1;
incrlt=1;
kplimit=3;

LONG = [-180:incrlat:180];
LT =   [6:incrlt:18];%LT bins 

% The reason for the negetive dip at the center of ACE stack for Kp < 2 is
% that since kp is low for the middle part,as comapred to earlier or later
% day/time, the center part will have low ampliutde of excursion than the
% edges. (vice versa for the center peak for kp > 6 (watsed 2 days to find
% out this !)




for ki = LONG(1):incrlat:LONG(end-1),
L = data(:,2) >= ki & data(:,2) <= ki+incrlat & ltt >LT(1) & ltt < LT(end) & data(:,10) <=kplimit ; %The last filter is to select only data before Dec 31, 2002

mydata = data(L,:);
mylt=ltt(L);
nlt=1;
for lt_cej=LT(1):incrlt:LT(end-1),
    LL =  mylt > lt_cej & mylt < lt_cej+incrlt & mydata(:,5) < 0;
%    cej_array(nbin,nlt) = sum(mydata(LL,5)<0);
    cej_array(nbin,nlt) = sum(LL);
       
    clt = mylt(LL);
    cut = mydata(LL,1);
    cdt = mydata(LL,5);

    
    numd = 1;
    cpr=0;

    for ii = 1:length(cut),
    lodp = find(data(:,1)==cut(ii));
    if (cut(ii)>cpr),
        if mean(diff(data(lodp-5:lodp+5,1))) > 0.06,
               cdata(nbin,nlt,numd,:) = data(lodp-5:lodp+5,5);
               ctime(nbin,nlt,numd,:) = data(lodp-5:lodp+5,1);
               numd=numd+1;
        else,
            cdata(nbin,nlt,numd,:) = data(lodp-10:2:lodp+10,5);
            ctime(nbin,nlt,numd,:) = data(lodp-10:2:lodp+10,1);
               numd=numd+1;
        end;
            cpr = data(lodp+5,1);
        end;
    end;
    
     
    
    nlt=nlt+1;
end



%fprintf('begin lat = %d end lat = %d number of cej = %d\n',ki,ki+incrlat,cej_data(nbin));
nbin=nbin+1;
end;

AA = [LONG(1:end-1)+diff(LONG)/2];

BB = [LT(1:end-1)+diff(LT)/2];
figure(1);
imagesc(AA, BB ,cej_array' )

 set(gca,'XTick',AA);
 set(gca,'XTickLabel',AA);
    
 
 figure(3);
 
 nn=1;
for i = 1:length(LONG)-1,
for j = 1:length(LT)-1,
subplot(12,12,nn);
plot(mean(squeeze(cdata(i,j,:,:))));
axis([-inf inf -inf inf]);nn=nn+1;axis off
end;
end;